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ABSTRACT 

The critical star formation rate density required to keep the intergalactic hydrogen 
ionised depends crucially on the average rate of recombinations in the intergalactic 
medium (IGM). This rate is proportional to the clumping factor C = {p'^)iGM/{Ph)^, 
where pb and (pb) are the local and cosmic mean baryon density, respectively and the 
brackets ()igm indicate spatial averaging over the recombining gas in the IGM. We 
perform a suite of cosmological smoothed particle hydrodynamics simulations that in- 
clude radiative cooling to calculate the volume- weighted clumping factor of the IGM at 
redshifts z ^ 6. We focus on the effect of photo-ionisation heating by a uniform ultra- 
violet background and find that photo-heating strongly reduces the clumping factor 
because the increased pressure support smoothes out small-scale density fluctuations. 
Photo-ionisation heating is often said to provide a negative feedback on the reionisa- 
tion of the IGM because it suppresses the cosmic star formation rate by boiling the 
gas out of low-mass halos. However, because of the reduction of the clumping factor 
it also makes it easier to keep the IGM ionised. Photo-heating therefore also provides 
a positive feedback which, while known to exist, has received much less attention. 
We demonstrate that this positive feedback is in fact very strong. Using conservative 
assumptions, we find that if the IGM was reheated aX z >9, the observed population 
of star-forming galaxies at z Ri 6 may be sufficient to keep the IGM ionised, provided 
that the fraction of ionising photons that escape the star-forming regions to ionise the 
IGM is larger than ~ 0.2. 

Key words: cosmology: theory - methods: numerical - hydrodynamics - radiative 
transfer - intergalactic medium - galaxies: formation 



1 INTRODUCTION 

The absence of a Gunn-Peterson trough in the majority 
of the observed absorption spectra towards high-redshift 
quasars suggests that the reionisation of intergalac- 
tic h ydrogen was complete d at a redshift z > 6 (see 
iFan. Carilli.fc Keating l2006l for a recent review) and that 
it remained highly ionised ever since. Current observational 
estimates of the ultr a-violet (UV) luminosity density at 
redshifts z < 6 (e.g. Stanway. Bunker, fc McMahon 20031: 



Lehnert fc Bremer 2003; Bunker ct al. 2004; Bouwcns ct al 
20041: lYan fc Windhorst 2004; Sawicki fc Thompson 2008: 



Stiavelli. Fall, fc Panagial |2004 iMalhotra erall |2005| ; 
Panagia et al.l 120051 ) . Taken at face value, these low SFR 
densities pose a severe challenge to commonly employed 
theoretical models in which the observed population of 
star-forming galaxies is the only source of ionising radiation 
in the high-redshift Universe. 

There are, however, large uncertainties associated with 
both the observ ationally inferred (see , e.g., the comprehen- 
sive analysis of iBouwens et al.l 120071 ) and the critical SFR 
densities. The critical SFR density. 



Bouwens et al. l2006l: [Mannucci et all |2007| ; lOesch et all 
20081 : iBouwens et al.l l2008r ). on the other hand, may 
imply star formation rate (SFR) densities several times 
lower than the critical SFR density required to keep 
the intergalactic medium (IGM) ionised (but see, e.g.. 
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here rescaled to match the most recent WM AP estimate for 
the cosmic baryon density |Komatsu et al.l 12008). has been 
derived by iMadau. Haardt. fc Rees (119991) using an early 
version of the lBruzual fc CharlotI 1 20031 ) population synthe- 
sis code, assuming a Salpeter initial stellar mass function 
(IMF) and solar metallicity. It results from simply equat- 
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ing the spatially averaged rate at which ionising photons 
are emitted into the IGM to the spatially averaged rate at 
which the intergalactic gas recombines. Eq. [T] is therefore 
incapable of addressing a number of potentially important 
physical effects. Some ionising photons will, for instance, 
redshift below the ionisation threshold before ionising and 
some ionising photons will have been emitted longer than a 
recombination time ago upon impact with a neutral atom, 
so that equating instantaneous rates is not appropriate and 
one even may have to take source evolution into account. It 
is therefore important to keep in mind that Eq. [l] is likely 
only accurate within factors of a few. 

The critical SFR is inversely proportional to the escape 
fraction fcsc, i.e. the fraction of ionising photons produced 
by star-forming galaxies that escape the interstellar medium 
(ISM) to ionise the IGM, and proportional to the average 
recombination rate in the IGM. The latter is parametrised 
using the dimensionless clumping factor C = {pb)iGM/(pb) , 
where pb is the baryon density, (pb) is the mean baryon 
density of the Universe and the brackets {)igm indicate 
spatial averaging over the gas constituting the recombining 
IGM. Under the assumption of a uniformly ionised IGM, the 
clumping factor expresses the spatially averaged number of 
recombinations occurring per unit time and unit volume in 
the ionised IGM, relative to that in gas at the cosmic mean 
density (pb)- A larger escape fraction implies a smaller crit- 
ical SFR density, as more photons are available to ionise the 
IGM. On the other hand, a larger clumping factor implies a 
larger critical SFR density since more ionising photons are 
required to compensate for the increased number of recom- 
binations. 

Most observational studies that comparqj the SFR 
density derived from estimates of the UV luminosity den- 
sity at redshift z ~ 6 to the critical SFR density as- 
sume an escape fraction /esc ^ 0.5 and a clumping fac- 
tor C = 30. While a v ariety of both observationa l and 



theoretical studies (e. g. Ilnoue. Iwata. fc Deharveu: 
and references therein; Razoumov fc Sommer-Larse: 



l2006l 
l2006l : 



iGnedin. Kravtsov. fc Chenll2008al ) have ruled out large 
cape fractions, the estimate for the clumping factor comes 
from a single cosmological simulation performed more than 
10 years ago IjGnedin fc Ostrikerlll997l ). It is on the basis of 
these values for the escape fraction and the clumping factor 
that the observed population of galaxies has been found to 
be incapable of keeping the intergalactic hydrogen ionised, 
forming massive stars at a rate which is up to an order of 
magnitude lower than required by Eq. [T] 

It has been pointed out that this discrepancy 
between the inferred and critical SFR densities 
could be resolved if the employed clumpin g fac- 
tor were too high (e.g. ISawicki fc ThompsonI 120061 : 
see also the discussion in iBouwens et al.l 120071 ). In- 
deed, in most (but not all) of th e more recent 
theoretical theoretical studies (e.g. IValageas fc Silkl 



^ Note that although the critical SFR density is sensitive to the 
IMF, this comparison is insensitive to the IMF provided the 
same IMF is used to compute the critical and observationally 
derived SFR densities. This is because the UV luminosity den- 
sity is dominated by the same massive stars that are responsi- 
ble for the emission of ionising photons with energies > 13.6 eV 
llMadau. Haardt. fc Ree3ll999fi . 
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clumping factors were derived. On the other hand, it is 
sometimes emphasised that simulations underestimate the 
clumping factor, due to a l ack of resolution (see, e.g., 
iMadau. Haardt. fc Reed Il999l ). In this work we perform 
a set of cosmological Smoothed Particle Hydrodynamics 
(SPH) simulations that include radiative cooling and photo- 
ionisation by a uniform UV background in the optically 
thin limit to study the clumping factor of the IGM. 

We focus on the effect of photo-ionisation heating on 
the evolution of the clumping factor. Previous investiga- 
tions of the impact of photo-heating on the reionisation 
of the IGM have almost exclusively come to the conclu- 
sion that it acts as to provide a negative feedback. Photo- 
heating boils the gas out of the potential wells of dark 
matt er (DM) h alos with virial temper atures Tyir < 10^ K 
fe.g. iThoul fc W einber g 1996; Navarro fc Steinmet j [l997l: 
' Barkana fc Loebl [l999l ; iKitavama fc Ikeuchil I2OO0I: iGnedinI 



2000bl: lDiikstraetal.1 [20041: IShapiro.Iliev.fc Ragal 12004 
Hoeft et al.1 I2OO6I: ICrain et all l2007l: iMesinger fc Diikstra 



2008'; 'O kamoto. Gao. fc TheunsI I2OO8I : iPawlik fc Schavej 
2008). This inhibits the formation of stars in these low-mass 
halos and thus decreases the ionising emissivity, which makes 
it more difficult to reionise the Universe. The same mecha- 
nism that reduces the number of ionising photons that are 
emitted into the IGM does, howe ver, also affect the evolution 
of the clumping factor (e.g. iHaiman. Abel, fc Madaul 2001 



or the clumping lactor (e.g. itlaiman. Abel, fe Madaul IzUU II: 
Oh fc Haimanll2003l: iKuhlen fc Madaul l2005l: IWise fc Abell 



Furlanetto. Oh, fc Briggj I2OO6I : ICiardi fc Salvaterral 



In this paper we demonstrate that photo-heating sig- 
nificantly lowers the clumping factor and hence the average 
recombination rate in the IGM. While photo-ionisation heat- 
ing undoubtedly impedes the production of ionising photons, 
our results imply that it also makes it much easier to keep 
the IGM ionised. 

The paper is structured as follows. In Section [2] we give 
a detailed description of our set of simulations. In Section [3] 
we use our simulations to compute the clumping factor of 
the IGM. Finally, in Section (4] we discuss our results and 
their implications for the value of the critical SFR density. 



2 SIMULATIONS 

We use a m odified version of the N-body/TreePM/SPH code 
GADGET-2 (|Springelll2005n to perform a suite of cosmological 
SPH simulations including radiative cooling. 

The initial particle positions and velocities are obtained 
from glass-like initial cond itions using CMBFAST (version 4.1; 
ISeliak fc Zaldarriagalll996l ') and employing the Zeldovich ap- 
proximation to linearly evolve the particles down to redshift 
z — 127. We assume a flat ACDM universe and employ the 
set of cosmological parameters [flm , ^b , ^a , o"8 , JT-s , ft] 
given by [0.258,0.0441,0.742,0.796,0.963,0.719], in 
agreement with th e WMAP 5-year observations 
iJKomatsu et al.l [2008J). For comparison, we also per- 
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Table 1. Simulation parameters: comoving size of the simulation box, Lbox (default value: 6.25 h~^ Mpc); number of DM particles, N^^ 
(default value: 256'^); mass of dark matter particles, m,(jm (default value: 8.6 X 10^ h~^ Mq); additional reheating energy per proton, Cr 
(default value: 2 eV); reheating redshift, z^ (default value: 9); kinetic feedback from supernova winds, winds (default: no); cosmological 
parameters, WMAP (default: 5-year). The number of SPH particles initially equals N^^ (it decreases during the simulation due to star 
formation). Bold font indicates our default simulation. 
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Figure 1. Thermal evolution of gas with overdensity A = 1 for 
characteristic choices of the reheating parameters Zr and e^ , which 
are listed in Table [T] for the simulations indicated in the legend. 
Note that even in the absence of an additional energy input at 
redshift z = Zr, i.e. for Cr = 0, as it is the case for simulation 
r9L6N256lowT, the gas is quickly heated by the UV background 
to a temperature T ~ 10* K. 



form some simulations employing the set of cosmological 
parameters [0.238,0.0418,0.762,0.74,0.951,0.73] and 
[0.25,0.045,0.75,0.9. 1, 0.73], consistent with WMAP 3-year 
ilSpergel et all l2007i ) and WMAP 1-year (|Spergel et al.1 
120031 ) observations, respectively. Data is generated at 50 
equally spaced redshifts between z = 2Q and z = Q. The 
parameters of the simulations employed for the present 
work are summarised in Table [T] 

The gravitational forces are softened over a length of 
1/25 of the mean dark matter int er-particle distance. We 
emplo y the star formation recipe of lSchave fc Dalla Vecchial 
l|2008h . to which we refer the reader for details. Briefly, gas 
with densities exceeding the critical density for the onset of 
the thermo-gravitational instability (hydrogen number den- 
sities TiH ~ 10~^ — 10~^ cm~3) is expected to be multiphase 



and star- forming (|Schavell2004l ). We therefore impose an ef- 
fective equation of state (EoS) with pressure P oc p"*""'* for 
densities nn > n^, where ny = 10~^ cm"^, normalised 
to P/k — 10'^ cm~3 K at the critical density n^j. We use 
7cff = 4/3 for which both the Jeans mass and the ratio of 
the Jeans length and the SPH kernel are independent of 
the density, thus preventing spurious fragmentation due to 
a lack of numerical resolution. Gas on the effective EoS is 
allowed to form stars using a pressure-dependen t rate that 
repro duces the observed Kennicutt-Schmidt law (|Kennicuttl 
Il998l ). renormalised by a factoqj of 1/1.65 to account for the 
fact that it assumes a Salpeter IMF whereas we are using a 
Chabrier IMF. 

The gas is of primordial composition, with a hydro- 
gen mass fraction X = 0.752 and a helium mass fraction 
Y = 1 — X. Radiative cooling and heating are included as- 
suming ionisation equilibrium, using tables generated with 
the publicly available package cloudy (ve rsion 05.07 of the 
code last described by Ferland et al.lll998r ). as described in 
IWiersma. Schaye fc SmithI ( 20081 ) . The gas is allowed to cool 
by coUisional ionisation and excitation, emission of free-free 
and recombination radiation and Compton cooling off the 
cosmic microwave background. 

We perform a set of simulations including photo- 
ionisation by a uniform UV background in the optically 
thin limit at redshifts below the reheating redshift z^. These 
simulations are denoted with the prefix r (see Table [TJ. 
To study the effect of reionisation reheating, we compare 
these simulations to a simulation that does not include 
photo-ionisation {L6N256). Note that the photo-ionisation 
changes the density of free electrons and the ionic abun- 
dances. Both the cooling and heating rates are t herefore af- 
fected by the inclusion of a UV backgr ound (e.g. lEfstathioul 
1 19921 : IWiersma. Schaye fc Smit"hll2008l ) . 



■^ This conversion factor b e tween SFRs has been computed us- 
ing the iBruzual fc Charloti 1 120031) population synthesis code for 
model galaxies of age > lO'^ yr forming stars at a constant rate 
and is insensitive to the assumed metallicity. 
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Figure 2. Slices (of thickness 1.25 h~^ comoving Mpc) through the centre of the simulation box, showing the SPH overdensity field 
in the simulations L6N256 and r9L6N256 at redshifts z = 9.08 (left-hand panel; where they are identical) and 2; = 6 (middle panel: 
L6N256, right-hand panel: r9L6N256). The inclusion of photo-heating in r9L6N256 leads to a strong smoothing of the density field 
(right-hand panel). 



The properties of the UV background depend on the 
redshift of reheating. If z^ ^ 9, we employ the evolving 
UV background f rom q uasars and galaxies tabulated by 
iHaardt fc Madaul (I2OOI ) for z ^ z,. If z, > 9, we use the 

9 iHaardt fc Madaul (|200ll ) UV back ground for all red 



shifts 9 < z ^ Zi, and employ the evolving IHaardt fc Madaul 
l|2001l ) UV ba ckground for redshifts z SC 9. This is necessary 
because IHaar dt & Madau (2001) onl y tabulate up to z = 9 . 
For z > Zr, we employ the z ^ 9 IHaardt fc Madaul (120011 ) 
UV background but with its intensity at energies equal to 
and larger than 13.6 eV set to zero. Molecular hydrogen and 
deuterium and their catalysts are kept photo-dissociated by 
this soft UV background at all redshifts and therefore never 
contribute to the cooling rate. Our approach is motivated in 
the context of reionisation because the weak UV background 
established by the very first ionising sources is already suf- 
ficient to effici ently suppress the formatio n of molecular 
hydroge n (e.g. piairnan. Rees. fc Loeblll997| and referen ces 
therein; IClovedboOTl ; IChuzhov. Kuhlen. fc Shapiroll2007l ) . 

The reheating redshift z^ is a parameter in our sim- 
ulations. The most recent determination of the Thomson 
optical depth towards reionisation from the WMAP (5-year) 
experiment implies a reionisation redshift Zreion = 11.0±1.4, 
assuming that the transition from t he neutral to the full y 
ionised Universe was instantaneous (iKomatsu et al.l 120081 ). 
The Thomson optical depth towards reionisation provides, 
however, only an integral constraint on the Epoch of 
Reionisation. The reionisation history may therefore have 
been considerably more intricate. An early population 
of X-ray sources, for example, could reheat the IGM 



to temperatures 
shifts (e. , 
19991; \a 



10 K al ready at much higher red- 



CoUin-Souffrinl Il99ll: iMadau fc Efstathioul 
200 ll; Venkatesan. Giroux 



2001 



2004 



Shui: 
Machacek, Bryan, fc Abe l 2003; Madau et alj 
Ricotti fc Ostrikerl [200J). We therefore study a range 
of thermal histories, performing simulations using 
Zr = 7.5, 9, 10.5, 12, 13.5, 15 and 19.5. To be conserva- 
tive, we use the relatively low reheating redshift Zr = 9 as 
our default value. 

In our simulations we compute the photo-heating rates 
in the optically thin limit, which means that we under- 



esti mate the temperature o f the IGM during reionisation 
(e.g. lAbel fc Haehnelti 119991 ). We therefore inject an addi- 
tional thermal energy £ r per proton at z = Zr (see, e.g., 
iThoul & Weinberg 1996). By varying the parameter tr, we 
will investigate the sensitivity of our results to the temper- 
ature of the reheated IGM. Our default simulation employs 
Er = 2 eV. Fig. [1] shows the thermal evolution of gas at 
the cosmic mean baryon density (pb), i.e. of gas with over- 
density A = Ph/{Pb) ~ 1, for different values of e^ and Zr- 
At z = Zr, the gas is heated to Tr ~ 10* K for e^ ~ '2 eV, 
whereas the gas temperature is about an order of magnitude 
higher (lower) for e^ = 20 eV (er = eV). After reheating 
the gas quickly looses memory of its initial temperature and 
by 2 = 6 the gas temperature is T ~ lO'* K in all cases. 

In one of our simulations {r9L6N256winds) we in- 
clude kinetic fe edback from star format i on. W e employ the 
prescription of Dalla Vecchia fc Schaye (^200»), which is a 
variation of the Springel fc Hernquist (2003) recipe for ki- 
netic feedback. In this prescription, core-collapse super- 
novae locally inject kinetic energy and kick gas particles 
into winds. The feedback is specified by two parameters, 
the mass loading 77 = M^/M,, which describes the ini- 
tial wind mass loading M^ in units of the cosmic SFR 
Mt, and the initial wind velocity v^n- We use rj = 2 and 
600 km s~^, consistent with observations of local 



(e.g. IVeilleux. Cecil fc Bland-Hawthornll2005l ) and redshift 
z ~ 3 ( Shaplev et al.l |2003| ) starburst galaxies. Note that 
wind particles are not hydrodynamically decoupled and that 
they are launched local to the star for mation event, different 
from the ISpringel fc Hernouisti l|2003l ) recipe. 



3 RESULTS 

In this section we employ the set of simulations described 
in Section [5] and summarised in Table [T] to calculate the 
clumping factor of the IGM. We start in Section 13.11 with 
analysing the distribution of the gas in our default simu- 
lation r9L6N256 and in the simulation L6N256, which is 
identical to our default simulation except for the fact that it 
does not include a photo-ionising background. In Section [3^ 
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we discuss the definition of the clumping factor and com- 
pare the clumping factors derived from our default simula- 
tion r9L6N256 to that derived from simulation L6N256. We 
discuss the convergence of our results with respect to varia- 
tions in the mass resolution and in the size of the simulation 
box in Section r3.2.1l In Section r3.2.2l we vary the redshift at 
which the ionising UV background is turned on and in Sec- 
tion 13.2.31 we demonstrate that our conclusions are robust 
with respect to our choice for the temperature to which the 
IGM is photo-heated. In Section r3.2.4l we discuss how kinetic 
feedback from supernova winds affects our results and quote 
the clumping factors obtained from the simulations employ- 
ing WMAP 3-year and 1-year cosmological parameters. We 
conclude with a brief comparison to previous work. 

3.1 The gas density distribution 

Here we compare the gas distributions in our default sim- 
ulation r9L6N256 (in which the UV background is turned 
on at redshift Zr = 9) and in the simulation L6N256 (which 
does not include photo- ionisation heating). 

Figure [2] shows the overdensities at redshifts z = 9.08 
and z = 6 in a slice through the simulation box for these 
simulations. Heating by the photo-ionising background al- 
most instantaneously increases the gas temperatures to 
Tr ~ 10* K (see Fig. [l| and accordingly raises the cos- 
mological Jeans mass. Gas that had already settled into 
the potential wells of DM halos with virial temperature 
jTvir ^ 2r is driven back into the diffuse IGM by the in- 
creased pressure gradient for ces (e.g. lBarkana fc Loeblll999l : 
IShapiro. Iliev. fc Ragall2004 ). The large cosmological Jeans 
mass prevents any re-accretion of gas or infall of previ- 
ously unbound gaseous material into these low-mass halos 
and keeps the IGM diffuse. Comparing the middle panel 
with the right-h and panel of Fig. (2] this J eans fi l tenni 
(e.g. [Sha piro. Gir oux. fc Ba bul 1994; Gncdin fc Huilll99 ' 
iGnedinllioOOb, : .Okamoto. Gao. fc Theuns„2008i ') in the pres 



ence of photo-heating leads to a strong smoothing of the 
small-scale density fluctuations by z = 6. 

A detailed analysis of the overdensity distribution in 
the simulations can be obtained by studying ■Pv(A), the 
volume-weighted probability density function (PDF) for A. 
We show the PDF (per unit logj^g A and normalised accord- 
ing to /y°° dA Pv(A) = 1) in the left-hand panel of Fig. [3 It 
is important to be aware of the fact that the finite numerical 
resolution of our simulations implies an unavoidable intrinsic 
smoothing of the gas density distribution on the scale of the 
SPH kernel or the scale over which the gravitational forces 
are softened, whichever is larger. A numerical smoothing on 
scales larger than the Jeans filtering scale (below which the 
gas density distribution is physically smoothed) would pre- 
vent us from obtaining converged results. We will discuss the 
convergence of our simulations with respect to resolution in 
Section 13.2.11 

At redshift z — 9.08 (black solid histogram), the gravi- 
tational amplification of the overdensities present in the ini- 
tial conditions has produced a significant deviation of the 
PDF from its primordial Gaussian shape. The flattening of 
the slope of the PDF for overdensities 1 < logj^g A < 2 can 
be attributed to the shock-heating of gas falling into the po- 
tential wells of dark matter halos, most of which have virial 
temperatures < 10* K, which we refer to as low-mass halos. 



The shape of the PDF is determined by the effective EoS 
once the gas reaches the critical density for the onset of star 
formation (njj = 10^^ cm~^, see Section[51 indicated by the 
vertical lines). 

At redshifts z < Zr, the shape of the PDF strongly de- 
pends on whether photo-heating by the ionising background 
is included or not. In the absence of such a background 
{L6N256, blue dashed histogram), gravitational collapse 
proceeds unimpeded, increasing the PDF at logj^Q A > 1. 
Since the gas that accretes onto DM halos must originate 
from the reservoir at logj^g A < 1 (the diffuse IGM), the 
PDF decreases over this range of overdensities. As a result, 
the maximum of the PDF shifts to lower overdensities. At 
the same time, the gravitational amplification of large-scale 
underdense regions leads to an increase in the PDF around 
overdensities logjQ A ~ — 1. 

Photo-heating in the presence of the ionising back- 
ground photo-evaporates the gas in DM halos, as described 
above. The bump in the PDF around 1 < logjQ A < 2 
therefore disappears (red dot-dashed histogram). Note that 
the redistribution of the baryons due to photo-heating also 
slightly increases the minimum overdensity that is present 
in the simulation. In Appendix IA2I we compare the PDF 
ob tained from our default simulation to th e fit provided 
by iMiralda-Escude. Haehnelt. fc Reed l|200Ci ). which is of- 
ten employed to compute the clumping factor in (semi- 
)analytical studies of the epoch of reionisation. 



3.2 The clumping factor 

In this section we demonstrate how the clumping factor, 
C = {pb)iGM/{pb)^, depends on the definition of the IGM 
and compute it for our default simulation r9L6N256 and for 
the simulation L6N256. This allows us to investigate how 
the clumping factor is affected by the inclusion of a photo- 
ionising background. 

Our main motivation for computing the clumping fac- 
tor of the IGM is to evaluate the critical SFR density re- 
quired to keep the IGM ionised. The critical SFR density 
describes the balance between the number of ionising pho- 
tons escaping into the IGM (parametrised by the escape 
fraction) and the number of ionising photons that are re- 
moved from the IGM due to photo-ionisations of recombin- 
ing hydrogen ions (parametrised by the clumping factor). 
When the ratio of photon escape rate to recombination rate 
is larger than unity, the rate at which galaxies form stars ex- 
ceeds the critical SFR density and is thus sufficient to keep 
the IGM ionised. 

It is important to realise that only recombinations lead- 
ing to the removal of ionising photons which escaped the 
ISM of the star-forming regions contribute to the balance 
that gives rise to the definition of the critical SFR den- 
sity. To separate the gas in the ISM from the gas in the 
IGM , a simple threshold density criterion is of ten employed 
(e.g. iMiralda-Escude. Haehnelt. fc Reesll2000l : see also the 
discussion in iMiralda-Escudg l2003h . Ionising photons are 
counted as escaped once they enter regions with gas densities 
pb < Pthr. Consequently, only gas with densities pb < Pthr, 
or equivalently, gas with overdensities A < Athr = pthr/(pb) 
should be considered in the evaluation of the clumping fac- 
tor. 

The threshold density pthr depends on which gas is con- 
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Figure 3. Left-hand "panel: Volume-weighted PDF of the baryon overdensity A per unit logjQ A for simulations r9L6N256 and L6N256 
at rcdshifts z = 9.08 and ^ = 6, as indicated in the legend. Photo-heating destroys the bump around overdensities 1 < logjQ A < 2, which 
mark the gas that accretes onto DM halos. The vertical lines (which match the colour and style of the corresponding PDFs) indicate 
the overdensities corresponding to the onset of star formation. Right-hand panel: Clumping factor C(< Atij^) of gas with overdensity 
A < Atijj. for the simulations shown in the left-hand panel. The inclusion of photo-heating in r9L6N256 leads to a clumping factor 
that is substantially smaller than that obtained from L6N256, for threshold overdensities An^^ > 10. Note that the maximum threshold 
overdensities we consider for the calculation of the clumping factor are given by the critical density n^ = 10^^ cm^'' for the onset of 
star-formation (the vertical lines shown in the left-hand panel). 



sidered to be part of the ISM, and which gas is considered 
to be part of the IGM. As long as the definition of the es- 
cape fraction and that of the clumping factor refer to the 
same decomposition of the gas into IGM and ISM, its value 
can be chosen arbitrarily. We therefore treat the threshold 
density as a parameter and compute the clumping factor as 
a fun ction of Athr (cp. iMiralda-Escude. Haehnelt. fc ReesI 
[2003), 



C(< At 



dA A^ Pv(A), 



(2) 



where 'Pv(A) is normalised according to J "" dA Pv(A) = 
1. In practice, we calculate C(< Athr) by performing a 
volume-weighted summation over all SPH particles with 
overdensities At < Athr, i-e 



C{< At: 



E 



A,<At 



ft A 



Ea, 



<Ath, 



hf 



(3) 



where hi is the radius of the SPH smoothing kernel associ- 
ated with SPH particle i. We verified that replacing /if with 
rrii/pi as an estimate for the volume occupied by SPH par- 
ticle i (with mass irn) gives nearly indistinguishable results. 
By definition, C(< Athr) increases monotonically with 
the threshold density pthr- Here, we set an upper limit 
to Athr, corresponding to the threshold density n^ = 
10~^ cm~^ for the onset of star formation that we employ 
in our simulations. Since we impose an effective EoS for gas 
with densities larger than rin (Section [2]), its PDF is not 
expected to reflect the PDF of real star-forming regions, 
motivating our choice for the maximum threshold density. 
The choice is conservative, leading to an overestimate rather 
than an underestimate of the critical SFR density, since the 
threshold density marking the escape of ionising photons 
and hence the clumping factor of the IGM to be used in 



Eq.[T]i s likely to be lower (see, e.g., the discussion in lGnedirJ 
l2008bl ). 



In the right-hand panel of Fig. Owe show C(< Athr) for 
the simulations r9L6N256 and L6N256 at redshifts shortly 
before (at z = 9.08, when r9L6N256 and L6N256 are iden- 
tical) and well after (at z = 6, when they differ by the pres- 
ence and absence of an ionising background, respectively) 
the reheating redshift z-r = 9. In agreement with our dis- 
cussion above, the clumping factor increases monotonically 
with the threshold density. Its dependence on redshift can 
be understood by looking at the evolution of the shape of 
the PDF, which we discussed in the previous section. 

For L6N256, i.e. in the absence of photo-heating, the 
clumping factor for threshold overdensities log^Q Athr > 1 
is larger at 2 = 6 than at 2; = 9.08, due mainly to the 
growth of the bump present in the PDF for overdensities 
1 < logj^Q A < 2. For logj^Q Athr ~ 0, on the other hand, the 
clumping factor is slightly smaller at 2 = 6 than at 2 = 9.08, 
which is caused by the depletion of the diffuse IGM through 
accretion onto DM halos. Note that at 2 = 6 the clumping 
factor reaches a maximum value of C « 20, which is signifi- 
cantly smaller than the value quoted bv lGnedin fc Ostriken 
(19971), which is commonly employed in observational stud- 
ies. This is probably because lGnedin fc Ostrikeij (J19971 ) com- 
puted the clumping factor including gas of any density, i.e. 
using a density threshold implicitly set by the maximum 
overdensity resolved in their simulation. 

The evolution of the clumping factor in r9L6N256, i.e. 
in the presence of the ionising background, is very different. 
At 2 = 6 it is close to that at z — 9.08 for all threshold over- 
densities. Compared to L6N256, the difference between the 
clumping factors for 2 = 6 and 2 = 9.08 is greatly reduced 
and the clumping factor at redshift 2 = 6 never reaches 
values larger than C ~ 6. 
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Figure 4. Clumping factor C{< Athr) of gas with overdensities 
A < Atiir and its dependence on resolution (at fixed box size, 
left-hand panel) and on box size (at fixed resolution, right-hand 
panel). The clumping factor obtained from our default simulation 
r9L6N256 is converged with respect to the employed resolution 
for all threshold overdensities shown. With respect to the size 
of the simulation box, it is converged for threshold overdensities 
log^Q Athr ^ 2. For larger threshold overdensities, full conver- 
gence may require the use of simulation boxes even larger than 
12.5 h~^ comoving Mpc, the size of the largest box employed 
here. 



3.2.1 Convergence tests 

In this section we check whether our results are converged. 
Generally, one expects the clumping factor to increase with 
both the spatial resolution and the size of the simulation 
box. The spatial resolution determines the smallest scale 
on which fluctuations in the density field may be identi- 
fied, whereas the size of the simulation box sets a cut-off 
to the largest scale on which the overdensity field can be 
non-zero. Moreover, the size of the simulation box limits the 
mass of the largest halo present in the simulation. Fig. |4] 
demonstrates that our default simulation (r9L6N256) is of 
sufficiently high resolution and employs a sufficiently large 
box to allow a faithful computation of the clumping factor of 
the reheated IGM at z = 6. In the left-hand panel we show 
the clumping factor in three simulations that use the same 
box size, but have mass resolutions that differ by multiples 
of 2'' , whereas in the right-hand panel we show the clumping 
factor in three simulations that employ the same resolution, 
but have box sizes that differ by multiples of 2. 

When compared to our default simulation r9L6N256, 
decreasing the mass resolution by factors of 8 {r9L6N128) 
and 64: (r9L6N064) only leads to insignificant and unsys- 
tematiqj changes in the clumping factor (left-hand panel). 
This can be understood by noting that the virial mass of 
halos corresponding to a virial temperature Tvir = Tr is re- 
solved with > 100 particles for redshifts z < 9. Any fur- 
ther increase in resolution would therefore mostly affect the 
abundance of DM halos with Tvir ^ Tr- In the presence 
of the UV background the gas in these halos is, however, 
quickly photo-evaporated (if formed at « > 2r) or prevented 
from accreting by the large Jeans mass associated with the 
photo- heated gas (for z < Zr). At redshifts well below the 



reheating redshift Zr, i.e. when sufficient time has passed to 
accomplish the photo-evaporation of halos and to allow the 
gas to respond hydrodynamically to the jump in the Jeans 
mass (i.e. to Jeans-filter the IGM), the observed convergence 
with respect to mass resolution is therefore expecteqj. 

Since the clumping factor has already converged for the 
spatial resolution used in simulation r9L6N128, we can em- 
ploy this simulation to verify whether the size of the box 
of our default simulation r9L6N256 is sufficiently large to 
enable an unbiased estimate of the clumping factor (right- 
hand panel). Overall, increasing the box size (at fixed res- 
olution) from 6.25 h~^ comoving Mpc by a factor of two 
to 12.5 h~^ comoving Mpc leaves the clumping factor al- 
most unaffected. For threshold densities log^ Athr ^ 2 the 
clumping factor may, however, not yet have fully converged, 
indicating that even larger simulation boxes than that con- 
sidered here ma y be required for its co mputation (see also 
the discussion in lBarkana fc Loebll2004 ). 



3.2.2 Varying the reheating redshift 

To study the effect of photo-heating on the evolution of the 
clumping factor in more detail, we make use of the clumping 
factors 



C_i = C(< 10"' cm-'mH/(X{/9b))) 



and 



Cioo = C(< min(100, 10"' cm-^mn/{X{p^)))). 



(4) 



(5) 



C_i is the clumping factor for gas with densities below 
riH ~ jih — 10^^ cm^^, the maximum threshold density 
we consider. For redshifts z < 16.3, Cioo fixes the thresh- 
old overdensity to a value between the me an overdensity 
of sp herical top-hat DM halos, (~ IStt^; e.g. IPadmanabharJ 
[1993) and the overd ensity at the virial r adius of an isother- 
mal DM halo (« 60; iLacev fc Colelll994 ) and closely agrees 
with the threshold densities commonly employed in the lit- 
erature to calculate the clumping factor of the IGM. For 
larger redshifts, the threshold overdensity Athr ~ 100 corre- 
sponds to a density tih > tih = lO"'^ cm""^, which is larger 
than the critical density for the onset of star formation in 
our simulations. The definition Eq.[S]ensures that the maxi- 
mum density that we consider for the computation of Cioo is 
always smaller than n^. Note that C_i refers to the clump- 
ing factor of gas with densities below a fixed proper density, 
while Cioo refers to the clumping factor of gas with densi- 
ties below a fixed overdensity for redshifts z < 16.3 and is 
identical to C_i for larger redshifts. 

The evolution of C_i and Cioo is shown, respectively, 
in the left-hand and right-hand panels of Fig. [S] for the sim- 
ulations r9L6N256 (i.e. reheating at redshift Zr = 9) and 
L6N256 (i.e. no reheating). In the same figure we also in- 
clude the evolution of C_i and Cioo obtained from the set 
of simulations r[zr]L6N256, where Zr = 7.5, 10.5, 12, 13.5, 15 
and 19.5. While the simulations are identical for 2 > Zr, the 
ionising background strongly affects the evolution of C_i(z) 



^ Note that the clumping factor is even slightly larger in 
r9L6N128 than in r9L6N256, although the latter has an 8 times 
higher mass resolution. 



■* We note that in the absence of a photo-ionising background, 
convergence may be more difficult to achieve, requiring a higher 
mass resolution than employed here. Convergence will be easier 
to achieve for lower threshold densities. 
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Figure 5. Evolution of tlie clumping factors C-i (left-hand panel) and Cioo (right-hand panel) for different reheating redshifts Zt, as 
indicated in the legends. Note that C—i(z) and Ciooiz) in all reheating simulations evolve towards the clumping factors obtained in 
the simulation with reheating at Zr = 19.5 (bottom black solid curves). At z = 6, the clumping factors are therefore insensitive to the 
reheating redshift provided that z^ ^ 9, with C— 1(2 = 6) ~ 6 and Cioo(z = 6) ~ 3. Fits to the evolution of the clumping factors are 
given in Appendix lAll 



and Cioo(2) for z < Zi. Whereas in L6N256 the clumping 
factors reach C_i « 20 and Cioo ~ 8 at z = 6, they only 
reach C_i ~ 6 and Cioo ~ 3 in r9L6N256, which are smaller 
than in L6N256 by roughly a factor of three. Note that be- 
cause the clumping factor obtained from simulation L6N256 
is likely to be not fully converged with respect to resolu- 
tion (see footnote Q, we may even have underestimated the 
magnitude of the decrease in the clumping factor due to 
photo-heating. 

Interestingly, the clumping factor at z = 6 is insensitive 
to the redshift z^ at which the UV background is turned on, 
as long as Zr > 9. This is because, after an initial transitory 
phase, the evolution of the clumping factor obtained for re- 
heating at redshift z,- approaches that obtained for reheating 
at Zt = 19.5. Note that the difference between the clump- 
ing factors obtained from r[zi]L6N256 and rl9.5L6N256 
becomes smaller with increasing reheating redshift z-^. In 
particular, the clumping factors obtained for reheating at 
Zr = 15 are nearly identical to those obtained for reheating 
at Zr = 19.5 at all redshifts. The evolution of the clumping 
factors obtained from rl9.5L6N256 can therefore be consid- 
ered to reflect the evolution of the clumping factor in the 
limit of reheating at very high redshift, z-^ 3> 19.5. 

In Appendix I All we provide fits to the evolution of the 
clumping factor over the redshift range 6 ^ z ^ 20 for a 
range of (over-)density thresholds and reheating redshifts. 
These fits (Eqs. lAll and IA2p may be employed in (semi- 
) analytic models of the epoch of reionisation. Many such 
models assume that reionisation heating provides only a neg- 
ative feedback on the reionisation process, reducing the star 
formation rate due to the photo-evaporation of gas in low- 
mass halos. However, as we have shown here, photo-heating 
decreases the clumping factor, and hence the average recom- 
bination rate. Since this makes it easier to keep the IGM 
ionised, reionisation heating also provides a positive feed- 
back on the process of reionisation. Although the relative 
importance of both can only be assessed using larger hydro- 



dynamical simulations of higher resolutiorjj, it is clear that 
models that do not account for this positive feedback will 
underestimate the efficiency with which star-forming galax- 
ies are able to reionise the IGM. 

3.2.3 Dependence on the reheating temperature 

In this section we investigate the robustness of our results 
with respect to our simplified treatment of photo-ionisation 
heating. 

We have considered photo-heating by a uniform ionising 
background in the optically thin limit. In reality, the reion- 
isation process is likely to be driven by inhomogeneously 
distributed sources in an initially optically thick medium. 
The temperature to which the IGM is reheated will then 
not only depend on the spectrum of the ionising sources, 
but also on the amount of spectral hardening due to the 
preferenti al absorption of the less energetic ionising pho- 
tons (e.g. lAbel fc Haehneltlll999l : iBolton. Meiksin. fc White! 
120041 ). Moreover, the speed at which a particular patch of 
the IGM is reionised determines the duration during which 
its gas can cool efficiently, as the cooling is dominated by in- 
elastic collisions between free electrons and neutral atoms. 
Different reionisation histories may therefo re result in differ- 
ent IGM temperatures (e.g. Miralda-Escude fc ReesI 1994 : 
Theuns et alll2002l : IHuI fc Haimanll2003l : iTittlev fc MeiksinI 
2007f ). 

To bracket possible scenarios, we have performed two 



^ At 2 = 6, the cosmic SFR density (the stellar mass) in 
r9L6N256 is smaller than that in L6N256 by a factor 1.26 (1.17). 
In our simulations, photo-heating thus decreases the SFR density 
less strongly than it reduces the clumping factor, which would 
imply that the positive feedback is more important. The SFR 
densities in both r9L6N256 and L6N256 are, however, not fully 
converged with respect to resolution and box size. A final as- 
sessment as to whether the negative or the positive feedback is 
stronger must therefore be deferred to future studies using simu- 
lations with higher resolution and larger box sizes. 
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Figure 6. The dependence of C—i (upper set of curves) and Cioo 
(lower set of curves) on the value for the reheating energy Cr • Both 
r9L6N256highT (e^ = 20 eV) and r9L6N256lowT (e^ = eV) 
give results very similar to that obtained from the default run, 
r9L6N256 (tr = 2 eV), which demonstrates the robustness of our 
conclusions with respect to changes in the reheating temperature. 



additional simulations in which we varied the amount of en- 
ergy transferred to the baryons during the photo-ionisation 
process, r9L6N256highT and r9L6N256lowT. Whereas in 
the former we employ an additional energy input that is 
ten times larger than our default value (tr ~ 20 eV), in the 
latter no additional energy is injected (er — eV). We show 
in Fig. |6]that the evolution of C_i and Cioo obtained from 
these two simulations is very similar to that obtained from 
our default run, r9L6N256. The dependence on the reheat- 
ing temperature Tr > 10* K is weak, because halos with 
virial temperatures Tvir < 10* K are already efhciently de- 
stroyed for Tr « 10* K. A further increase in the reheating 
temperature mostly affects the fraction of mass in halos with 
larger virial temperatures, which is small. Moreover, Fig. [T] 
shows that the gas in the simulations r9L6N256highT and 
r9L6N256lowT quickly loses memory of its thermal state at 
some higher redshift, which is another reason for the simi- 
larity in the results obtained using different values for the 
reheating energy tr. 



3.2.4 Effect of kinetic supernova feedback and dependence 
on cosmological parameters 

The inclusion of kinetic feedback from supernovae in 
r9L6N256winds only weakly affects the evolution of the 
clumping factors. At redshift z — 6, C_i and Cioo 
are slightly larger (by factors 1.1 and 1.18, resp.) in 
r9L6N256winds than in our default simulation, r9L6N256, 
which does not include kinetic feedback. The reason for the 
slight increase in the clumping factors is that winds move 
gas from regions of densities larger than the critical density 
for the onset of star formation to regions of lower density 
that contribute to the calculation of the clumping factors. 
We note that the inclusion of kinetic feedback does, on the 
other hand, strongly affect the cosmic SFR. At 2: = 6, the 
cosmic SFR (the stellar mass) is lower in the simulation that 
includes kinetic feedback {r9L6N256winds) than in our de- 
fault simulation {r9L6N256) by a factor of 6.1 (3.4). 



Finally, we quote the clumping factors obtained from 
the simulations r9L6N256W3 and r9L6N256Wl, which 
employed cosmological parameters consistent with the 
WMAP 3-year and 1-year observations, respectively. We find 
that at redshift z — 6, the clumping factors C_i and Cioo 
are larger in r9L6A^256 than in r9L6A''256W3 by factors of 
1.31 and 1.16, respectively. They are smaller in r9L6N256 
than in r9L6N256Wl by factors of 0.74 and 0.84. In sum- 
mary, with respect to r9L6N256, the clumping factors are 
larger in r9L6N256Wl and smaller in r9L6iV256W3, as ex- 
pected from the corresponding values of as, which set the 
average absolute amplitude of the overdensity fluctuations. 

3.2.5 Comparison with previous work 

We conclude our study of the clumping factor with a brief 
comparison with previous work, shown in Fig. [T] The evolu- 
tion of the clumping factors in our simulations L6N256 and 
r9L6N256 is shown by the black solid and red dashed curves, 
respectively, where the upper (lower) set of curves shows C-i 
(Cioo)- We comp are it to the evolution of the clumping fac - 
tor presented i n iMiralda-Escude. Haehnelt. fc ReesI (120001 ) 
and llliev et al.l (|2007l ). which are amongst the most com- 
monly employed works on the clumping factor and make 
use of sufficiently different techniques to bracket a range of 
possible cases. We caution the reader that such a direct com- 
parison is difficult and of limited validity because of the very 
diff erent assumptions underlying the in divid ual w orks . 

iMiralda-Escude. Haehnelt. fc liee3 (|2000l ') used 

the LIO hydrodyn a mical simulation presented in 
iMiralda-Escude et all (119961 ) to obtain the PDF of the 
gas density at redshifts z = 2, 3 and 4. The simulation 
was performe d using th e TVD hydrodynamical scheme 
described in I Ryu et al.l (119931 ). It used a box of size 
10 h~^ comoving Mpc, 144^ dark matter particles and 288'^ 
gas cells and employed cosmological parameters consistent 
with the first-year COBE normalization. The simulation 
included photo-heating from a uniform UV background, 
computed from the emissiviti es of the sources in the simu - 
lation. We refer the reader to iMiralda-Escude et al, 
for more details. iMiralda-Escude. Haehnelt. fc ReesI ( 



1996) 



2OOOD 



also provided fits to the gas density PDF and presented 
a prescription for its extrapolation to redshifts 2: > 4. We 
employed this prescription to compute the clumping factor 
evolution using Eq. [S] 

The evolutio n of the clumping factors C-i and C i qq ob - 
tained from the iMiralda-Escude. Haehnelt. fc ReesI (120001 ) 
PDFs is shown, respectively, by the top and bottom blue 
long-dashed curves. For redshifts 2; > 9, it closely agrees 
with the corresponding evolution obtained from our sim- 
ulation r9L6N256. For lower redshifts the agreement is 
less good, although the clumping factors never differ by 
more than factors ~ 2. The differences between their and 
our results are probably due to the use of different hy- 
drodynamical schemes, different cosmological parameters 
and different prescriptions for the UV background. The 
change in the slope of the clumping factor growth that 
can be seen at redhift z ~ 9 is likely due to the inclu- 
sion of photo-heating. That this change is much less pro- 
nounced than in our simulation r9L6N256 may be due to 
a more gradual build-up of the ionising background in the 
iMiralda-Escude. Haehnelt. fc ReesI (|2000l ) simulation. 
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llliev et alj (|2007l ) computed the clumping factor from 
a pure dark matter simulation. The simulation employed 
a box of size 3.5 h~ comoving Mpc and 1624 particles 
and was initialized with cosmological parameters consistent 
with the WMAP 3-year result^. The clumping factor was 
computed by averaging over all dark matter densities, and 
hence only implicitly makes use of an overdensity threshold 
(determined by the maximum overdensity present in their 
si mulation). We ref er the reader to the original description 
in llliev et al.l 1120071) f or more details. 

Iliev et al.l (|2007l ) provided the following fit to the evo- 



lution of the clumping factor in their simulation, 
Ciiiov07(z) = 26.2917exp(-Q.1822z-f 0.003505^^), 



(6) 



which is valid over the range 6 < 2 < 30. It is shown by the 
green dotted curve. Since it was derived from a pure dark 
matter simulation, Eq. |6]does not capture the hydrodynam- 
ical response due to reionisation heating. It should therefore 
be compared to the evolution of the clumping factor ob- 
tained from our simulation L6N256, which did not include 
photo-heating. A direct interpretation of such a comparison 
is, however, difficult, because Eq. |6]does not explicitly refer 
to an overdensity threshold. 

Our comparison clearly illustrates that there is a con- 
siderable spread in the clumping factor values quoted in the 
literature. The interpretation of many studies is complicated 
by the fact that they do not refer to a density threshold, 
which means that the result is determined by the numerical 
resolution of their simulations. 




Figure 7. Clumping factor evolution: comparison with previ- 
ous work. The black solid and red dashed curves are the clump- 
ing factors C_i (upper set of curves) and Cioo (lower set of 
curves) obtained from our simulations L6N256 and r9L6N256. 
The blue dashed curves show the evolution of the clumping fac- 
tors C_i (upper curve) and Ciqq (lower curve) derived from the 
[ :as de nsity PDFs presented in lMiralda-Escude. Haehnelt. fc Reea 
2000J) , which implicitly incorporate the effects of photo-heating. 
The green dotted curve shows the evolution of the clumping factor 
(defined without explicitly referring to a (over-)density threshold; 
instead, the overdensity threshold was set by the numer ical res- 
olution) the dark matter simulation of llliev et al.l | |2007T| . which 
does not include the effects of photo-heating. In both cases a di- 
rect comparison is difficult, because of the different assumptions 
underlying the individual works. 



4 DISCUSSION 

Several observational studies have claimed that the star for- 
mation rate (SFR) density at redshift 2 ~ 6 is smaller than 
the critical SFR density required to keep the intergalactic 
medium (IGM) ionised. In the absence of a large popula- 
tion of unseen sources of ionising radiation, this discrepancy 
between the two SFR densities would be in direct conflict 
with the high degree of ionisation inferred from the non- 
detection of a Gunn- Peterson trough in the majority of the 
line-of-sight spectra towards z <& quasars. 

The critical SFR density is inversely proportional to the 
spatially averaged fraction of ionising photons that escape 
into the IGM per unit time and proportional to the clumping 
factor C = (pb)iGM/{pb)^, a measure for the average recom- 
bination rate in the IGM. One may therefore ask whether 
the discrepancy between the observed and critical SFR den- 
sities could be resolved by changing the assumptions about 
the values of either of these two quantities. In this work we 
considered the hypothesis that most observational studies 
overestimate the critical SFR density because they employ 
a clumping factor that is too large. 

We re-evaluated the clumping factor, analysing the 
gas density distributions in a set of cosmological smoothed 
particle hydrodynamics simulations that include radia- 
tive cooling and photo-ionisation by a uniform UV back- 
ground in the optically thin limit. The clumping fac- 



° We note that they also computed the clumping factor in a sim- 
ilar simulation that was initialized with cosmological parameters 
consistent with the WMAP 1-year results. 



tor of the IGM depends critically on the definition of 
which gas is considered to bo part of the IGM. Following 
[Miralda-Escudc, Haehnelt, & Roes (200(J), we assumed that 
all gas with densities below a threshold density constitutes 
the IGM and computed the clumping factor as a function of 
this threshold density. In addition, we introduced two physi- 
cally well-motivated definitions, Cioo, the clumping factor of 
gas with overdensities A < 100 and C_i, the clumping factor 
of gas with proper densities below wh = n^ = 10~^ cm~"^, 
our threshold density for the onset of star formation. 

By comparing simulations that include photo-ionisation 
by a uniform UV background to one that does not, we 
showed that photo-heating strongly infiuences the evolution 
of the clumping factor of the IGM. Photo-ionisation heat- 
ing expels the gas from within halos of virial temperatures 
Tvir J; 10* K and prevents its further accretion by raising 
the Jeans mass in the IGM. By suppressing the formation of 
stars in these low-mass halos, photo-heating from reionisa- 
tion decreases the rate at which ionising photons are emit- 
ted into the IGM and is therefore correctly said to exert 
a negative feedback on the reionisation process. The fact 
that photo-heating also leads to a decrease in the clump- 
ing factor and hence provides a strong positive feedback by 
making it easier to keep t he IGM ionised, is however of- 
ten ov erlooked (but s ee, e.g.. 



20071) 



Haiman, Abel, & Madau 20011: 



Oh fc Haiman 200i iKuhlen fc Mada u 2005; Wise fc Ab3 



200£ ; iFurlanetto. Oh. fc Brig3 I2OO6I : ICiardi fc Salvaterri 



At redshift 2: = 6, we find that C_ 
and that these values are insensitive 



to 



6 and Cioo ~ 3 
the redshift 2,- 
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at which the UV background is turned on, as long as 
-2r ^ 9. These values for C_i and Cioo are at least three 
times smaller than they would be in the absence of photo- 
heating. We demonstrated that our default simulation is 
converged at z = 6 with respect to the employed res- 
olution. It is converged with respect to changes in the 
box size for threshold overdensities logj^Q Athr ^ 2. In Ap- 
pendix |X] we provide fits to the evolution of the clumping 
factor for various (over-)density thresholds and reheating 
redshifts. There we also compare the probability density 
function (PDF) of the gas densities at 2; = 6 obtained from 
our default simulation to the wid e ly use d fit provided by 
iMiralda-Escude, Haehnelt. fc Re3 (|2000l '). We update their 
fitting parameters to best fit the PDF from our default sim- 
ulation. Finally, we compared our results for the clumping 
factor of the IGM to those obtained in previous works. 

Since even our most conservative estimate for the 
clumping factor (C_i ~ 6) is five times smaller than 
the clumping factor that is usually employed to determine 
the capacity of star-forming galaxies to keep the z = 6 
IGM ionised, our results may have important implications 
for the understanding of the reionisation process. Setting 
C = 6 in Eq. [l] the critical SFR density becomes p, = 
0.005 fZc M0 yr~-^ Mpc~^. This is smaller than recent 
observational estimates for the SFR density at 2; ~ 6, 
p. = 0.022 ± 0.004 Mq yr"^ Mpc"^ (integrated to the 
observed 2 ~ 6 faint-end limi t L > 0.04 I/^^s and dust- 
corrected; |Bouwens_et_alJ|2007|), provided that /esc > 0.2. 

Our study thus suggests that the observed popula- 
tion of star-forming galaxies may be capable of keeping 
the IGM ionised, relaxing the tension between observation- 
ally inferred and critical SFR density in view of the ob- 
servation of a highly ionised IGM at redshifts z < 6. We 
note that a,t z ~ 7, the SFR density is estimated to be 
p* = 0.004 ± 0.002 M© yr"^ Mpc~^ (integrated to the 
observed z fa 7 faint-end limi t L > 0.2 ^2=3 and dust- 
corrected; |Bouwens_et_alJ[2008|), whereas the critical SFR 
density is p, = 0.008 f~c Mg yr"^ Mpc"^ (using C = 6). 
The observed population of star-forming galaxies at 2 « 7 
is therefore not able to keep the IGM ionised. If the Uni- 
verse were ionised by this redshift, then the sources that 
were responsible remain to be discovered. 

We caution the reader that the comparison of the crit- 
ical and observed SFRs is subject to considerable uncer- 
tainty. First, the SFR inferred from UV galaxy counts prob- 
ably underestimates the true SFR, because these counts miss 
UV galaxies fainter than the faint-end limit implied by their 
sensitivities. These galaxies may, however, significantly con- 
tribute to the UV luminosity density if the faint-end slope 
of the UV luminosity function is sufficiently steep. Comple- 
mentary estimates of the high-redshift star formation rate 
based on measurements of the high-redshi ft (2 = 4 — 7) 
gamma ray burst rate (JYiiksel et al.l I2OO8I ) and measure- 
ments of the Lyman-alpha fo rest opacity at redshifts 2 ~ 3 
l|Faucher-Giguere et al.l |2003) indeed suggest 2; ~ 6 SFRs 
that exceed those inferred from UV galaxy counts by fac- 
tors of a few. 

Second, the expression for the critical SFR (Eq. [1} is 
only approximate. As already mentioned in the introduction, 
this expression is based on equating the rate at which ionis- 
ing photons escape into the intergalactic medium to the rate 
at which the intergalactic gas is recombining, both averaged 



in space. It neglects effects like the cosmological redshifting 
of photons below the ionisation threshold energy and evo- 
lution of the ionising sources during a recombination time. 



Because the cross-section uhi ' 



for absorption of ionis- 



ing photons by neutral hydrogen decreases with increasing 
photon frequency v, these effects may become important for 
photons whose mean free path is comparable to the cosmic 
horizon. Our limit on the escape fraction required to keep 
the Universe at 2 ~ 6 ionised may therefore only be accurate 
within a factor of a few. 

Our simulations demonstrate that radiation- 
hydrodynamical feedback due to photo-ionisation heating 
plays a key role in shaping the properties of the IGM at red- 
shifts z > 6. We have studied the impact of photo-ionisation 
heating on the clumping factor of the IGM assuming a 
uniform ionising UV background in the optically thin 
limit. In reality the reionisation process will, however, be 
more complex. We have demonstrated the robustness of our 
conclusions with respect to uncertainties in the temperature 
of the IGM resulting from our simplified treatment of the 
reioinsation heating, but there are other factors whose 
importance is more difficult to assess. 

Our use of the optically thin approximation neglects 
self-shielding, a radiative transfer effect due to which ha- 
los that would otherwise be co mpletely photo-evaporated 
could keep some of their gas (e.g., Kitayama fc Ikeuchilbood : 
ISusa fc Umemurall2004 iDiikstra et al.ll2004D . Since the self- 
shielded gas remains neutral, it should be excluded when 
computing the clumping factor. Self-shielding becomes im- 
portant for A^Hi ^ 10^* cm~^, whic h for self-gra vitating 



gas clouds corresponds to densities (jSchavd I2OOII ) nn 



10" 



*(r/io- 



> 



where F is the HI photo- 



ionisation rate. In our discussion of the critical SFR we have 
conservatively adopted the clumping factor C_i, defined us- 
ing a threshold density nu = 10~^ cm~'^, which is an order of 
magnitude larger than the critical density for self-shielding. 

Self-shielding might affect our results because it may 
lower the sp eed with which halos a re pho to-evaporated. 
The work by llliev. Shap iro. & Raga (2005b) (in combina- 
tion with the work by LShapi ro, Ilicv, fc R aga 2004 ) shows 
that photo-evaporation times obtained in simulations that 
employ an optically thin UV background may differ by fac- 
tors of a few from those obtained in detailed radiation- 
hydrodynamical simulations. However, if any halos would 
resist photo-evaporation much longer than predicted by our 
approximate treatment of photo-heating, then the clumping 
factor of the IGM would be even lower, because self-shielding 
locks the gas that would otherwise contribute to the clump- 
ing factor of the IGM in its neutral state. 

If absorption by optically thick (self-shielded) clouds 
becomes important, then the mean free path of ionising 
photons may be set by the m ean distance between these 
clouds (IZuo fc PhinnevI 1 19931 ) instead of by the opacity 
of the diffuse IGM. In this case it might be appropri- 
ate t o supplement the iMiralda-Escude. Haehnelt. fc ReesI 
I2OO0I model for the computation of the average recombi- 
nation rate with a more di rect account for these clouds 
as discrete photon sinks ( e.g., llliev. Scannapieco. fc Shapird 
l2005al : ICiardi et aj] 12009 ). The effect of absorption by op- 
tically thick clouds is, however, largely degenerate with the 
ionising efficien cy of the populat ion of star-forming galaxies, 
as explained in llliev et al.l l|2007l ). If, on the one hand, these 
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clouds are ionising sources themselves, then their contribu- 
tion to the average recombination rate can be described by 
their escape fractions. If, on the other hand, these clouds do 
not host ionising sources, then their contribution to the aver- 
age recombination can be accommodated by changes in the 
ionising efficiency of star-forming galaxies because of their 
biased clustering around these galaxies. The effects of opti- 
cally thick clouds can therefore approximately be accounted 
for by adjusting the properties of the ionising sources. 

Self-shielding is only one example of the physical ef- 
fects that we are ignoring. The inclusion of metals and 
molecules, for instance, would increase the ability of the gas 
to c ool, which may lea d to an increase in the clumping factor 
(e.g. lMaio et al.ll2007h . In the presence of photo-heating and 
for the threshold densities we employ to compute the clump- 
ing factors we however expect this effect to be very small. A 
more efficient cooling due to metals and molecules may also 
enable star formation and associated kinetic feedback in low- 
mass haloes with virial temperatures Tvir < 10* K. Both are 
processes that we have ignored but which are likely to affect 
the clumping factor evolution fe.g. IWise fc Abelll2008l ). 

The evolution of the clumping factor will also depend on 
the morphology of the reionisation transition. Our treatment 
implicitly assumed that all gas with overdensities smaller 
than a threshold overdensity is uniformly ionised, while 
all gas with larger overdensities is fully neutral. This pic- 
ture probably only applies to the late stages of reionisa- 
tion, when individual ionised regions start to overlap and 
the only neutral gas that remains to be ionised is locked 
up in regions of high gas overdensities. Before overlap, 
other reionisation models may be more useful for the de- 
scription of the c lumpi ng factor evolution. For example, 
iFurlanetto fc OhI (|2005|) point out that if the large-scale 
dense regions are ionised first, the clumping factor may be 
somewhat larger than one would otherwise expect, because 
the photo ns are initially confined to these dense regions. 
Moreover, IFurlanetto. Haiman. fc Qhl (|2008l ') show that in 
fossil ionised regions, that is, regions in which the gas freely 
recombines, the clumping factor will generally be smaller 
than for regions in photo-ionisation equilibrium, because the 
densest gas which contributes most to the clumping factor 
becomes neutral first. 

The clumping factor is an important ingredient of (semi- 
)analytic treatments of reionisation. It would therefore be 
highly desirable to evaluate the approximations we have em- 
ployed in our simplified treatment of the photo-heating pro- 
cess using large high-resolution radiation-hydrodynamical 
simulations of reionisation that include cooling by metals 
and molecules and feedback from star formation. At the mo- 
ment such simulations are, however, not yet feasible. 

In fact, current state-of-the-art radiative transfer simu- 
lations typically make use of clumping factors in their "sub- 
grid" modules because they lack the resolution to resolve the 
dumpiness of the gas directly. In addition, they typically do 
not include the effect of photo-heating. In fact, many radia- 
tive transfer simulations ignore hydrodynamics altogether 
and assume the gas to trace the dark matter. 

With this work we hope to have presented a con- 
servative assessment of the clumping factor of the post- 
reionisation IGM that may provide a useful input to future 
(semi-)analytic models and simulations of the reionisation 
process. 
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APPENDIX A: FITTING FORMULAS 

Here we provide fits to the evolution of the clumping fac- 
tor for a range of overdensity thresholds and all reheat- 
ing redshifts considered for use with (semi-)analytical mod- 
els of reionisation. We also compare the probability den- 
sity function of the gas densities at redshift z — 6 ob- 
tained from our default simulation to th e fit provided by 
iMiralda-Escude. Haehnelt. fc Ree^ (|2000l '). 



Al Clumping factor 

In this section we fit the evolution of the clumping factors 
C-i and Cioo over the redshift range 6 ^ z ^ 20, based on 
the data presented in Fig. [S] In addition, we give fits to the 
evolution of the clumping factors C_2 and Ciooo, C500, C200 
and C50, where C-2 = C{< 10"^ cm"^mH/(X(pb))) and 
Ciooo = C{< min(1000, 10"^ cm-3mH/(X(pb)))) and sim- 
ilar for C500 , C200 and C50 . We first give fits to the evo- 
lution of the clumping factors for the simulations L6N256 
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Table Al. Parameter values a, (i, 7 and <5 to be used in Eg. lAll 
in order to fit the evolution of the clumping factors C—i and C-2 
obtained from the simulations L6N256 and rl9.5L6N256 (see 
also Fig. O. 





L6N256 


rl9.5L6N256 


a-i 


1.29 


1.21 


/3-1 


0.00 


-3.66 


7-1 


0.47 


0.00 


5-1 


5.76 


8.25 


a-2 


1.29 


1.16 


0-2 


0.00 


-2.47 


7-2 


0.44 


0.00 


<5-2 


4.68 


5.16 



and rl9.5L6N256, i.e. the simulations without reheating 
and with reheating at the highest redshift we considered. 
These fits are then used to obtain fits to the evolution of the 
clumping factors for reheating at the intermediate redshifts 
Zi — 7.5, 9.0, 10.5, 12, 13.5 and 15.0 by interpolation. 

We approximate the evolution of the clumping factors 
for the simulations L6N256 and rl9.5L6N256 by 



C{z) = 



-72+(5 



+ a, 



(Al) 



where C is either C_i, C_2, Ciooo, C500, C200, Cioo or C50 
and similar for a, /3, 7, 5. The values for the parame- 
ters a, /3, 7, 5 are Usted in Tables |Xl] and |M] The fit 
(Eq. IA1[I is accurate to within < 10%. We emphasize that 
it is only strictly valid over the fitting range 6 ^ 2; ^ 20. 
For Ciooo,C500,C2oo,Cioo and C50, i-e. for the clumping 
factors that are defined using an overdensity threshold, we 
forced, however, the fits to approach the correct high-z 
limit, i.e. C — > 1, by fixing a = 1 during the fitting pro- 
cedure. Note that the threshold densities used with C-i 
and C-2 correspond to threshold overdensities Athr < 1 
for redshifts z > 79.4 and z > 36.3, respectively and that 
Ciooo , C'soo , (^200 , Cioo and C50 become identical to C_i for 
redshifts z > 7.0,9.1, 12.7, 16.3 and 20.8, respectively. 

The values of C for reheating at redshifts Zi — 
7.5,9.0,10.5,12,13.5 and 15.0 (hereafter C') are fitted by 
interpolating between the fits (Eq. IA1[I to the evolution of 
the clumping factors obtained from L6N256 (hereafter C") 
and rl9.5L6N256 (hereafter C^^'^). That is, we write 



C"'(2) = W{Z)C\Z) + [1 - W{Z)]C^''-\Z), 



where 



w(z) = — 
^ ' 2 



erf 



-c 



+ 1 



and erf is the error function, 
2 



erf (2:) 



dz exp I 



(A2) 



(A3) 



(A4) 



The constants ("^' and r^' are listed in Tables IA3I and IA4I 
Eq.EHfits the data to within < 10%. 



A2 Probability density function 

iMiralda-Escude. Haehnelt. fc Reed (|2000l ) provided a con- 
venient four-parameter fit to the volume- weighted probabil- 



Table A2. Parameter values a, /3, 7 and 6 to be used in 
Eq. lAll in order to fit the evolution of the clumping factors 
ClOOOiC'soo, C'200, C'loo and C50 obtained from the simulations 
L6N256 and rl9.5L6N256 (see also Fig.|5}. 





L6N256 


rl9.5L6N256 


«1000 


1.00 


1.00 


/3iooo 


-1.00 


-2.89 


71000 


0.28 


0.00 


'5iooo 


6.29 


6.76 


0500 


1.00 


1.00 


/^SOO 


0.00 


-2.44 


7500 


0.34 


0.00 


(5500 


4.60 


5.68 


0200 


1.00 


1.00 


/3200 


0.00 


-1.99 


7200 


0.30 


0.00 


5200 


4.04 


4.49 


OlOO 


1.00 


1.00 


/3ioo 


0.00 


-1.71 


7100 


0.28 


0.00 


<5ioo 


3.59 


3.76 


"50 


1.00 


1.00 


/350 


0.00 


-1.47 


750 


0.23 


0.00 


550 


2.92 


3.08 



Table A3. Parameter values C,^' and r^' to be used in Ea. lA2l in 
order to fit the evolution of the clumping factors C— 1 and C-2 
obtained from the simulations r[z^]L6N256 (see also Fig. [Sj- 



Zy 


C\ 


rl\ 


Cl'2 


^2 


7.5 


6.83 


0.83 


6.61 


0.80 


9.0 


8.10 


1.25 


7.78 


1.26 


10.5 


9.41 


1.71 


8.92 


1.65 


12.0 


10.77 


2.16 


10.02 


1.96 


13.5 


12.10 


2.45 


11.09 


2.24 


15.0 


13.27 


2.43 


12.20 


2.48 



ity density function (PDF) of the gas density at redshifts 
2 = 2, 3 and 4 obtai ned from the LIO hydrodyna mical sim- 
ulation described in iMiralda-Escude et al.l (|l996l ) . In addi- 
tion, they provided a prescription for extrapolating this fit 
to higher redshifts. Here we compare the volume-weighted 
PDF obtained from our default simulation {r9L6N256) at 
2 = 6 to their fit. Because we find that their set of fitting pa- 
rameters provides a somewhat poor description of this PDF, 
we also provide an updated set of parameters that yields a fit 
which more accurately describes the one-point distribution 
of gas densities at 2 = 6 in our default simulation. 

In Fig. lAll we show the volume-weighted PDF of the 
gas density per unit logjQ A. We have investigated the con- 
vergence of this PDF with respect to the resolution and the 
size of the simulation box, using the same set of simulations 
that we employed to study the convergence of the clumping 
factor in Section [3.2.11 We find that the PDF is converged 
with respect to changes in the resolution. It is converged 
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Table A4. Parameter values (^^' and t^' to be used in Ea. lA2l in order to fit the evolution of the clumping factors CiooOi C'soOi ^2001 C'lOO 
and C50 obtained from the simulations r[zr]L6N256 (see also Fig. [5)l. 



'>50C 



*Zr 

''200 



AZr 

SlOO 



>z, 
'»50 



7.5 


6.71 


0.68 


6.74 


0.76 


6.71 


0.83 


6.60 


0.75 


6.53 


0.85 


9.0 


7.98 


1.29 


7.96 


1.22 


7.93 


1.32 


7.75 


1.72 


7.67 


1.37 


10.5 


9.49 


2.00 


9.50 


2.12 


9.33 


2.03 


8.99 


1.71 


8.94 


2.03 


12.0 


11.10 


2.58 


11.33 


3.05 


10.98 


2.90 


10.39 


2.38 


10.40 


2.84 


13.5 


12.36 


2.35 


13.12 


3.63 


12.75 


3.72 


11.93 


3.07 


12.02 


3.69 


15.0 


13.27 


1.77 


14.36 


3.26 


14.64 


4.49 


13.66 


3.75 


13.86 


4.60 



with respect to changes in the size of the simulation box for 
overdensities logjo ^ ^ 2. For larger overdensities, the val- 
ues of the PDF obtained from our default simulation may 
not yet be fully converged. We note that at z = 6 the PDFs 
obtained from the simulations r[zr]L6N256z that employ a 
reheating redshift Zr > 9 are almost undistinguishable from 
that obtained from our default simulation. 

We compar e the PDF to th e fit given by 
iMiralda- Escudc. Haehneh. fc Reed (|2000l 'l. 



Pv(A)dA = ylexp 



(A 



-2/3 



C,f 



A'^dA, 



(A5) 



2(25o/3)2 

using their values for the parameter^ A = 0.375, So = 
1.09, /3 = 2.50 and Co = 0.880. As can be seen from 
Fig. lAll for overdensities 1 < logj^Q A < 2 the PDF obtained 
from our default simulation is significantly steeper than pre- 
dicted by the fit. Using the simulations r9L6N256W3 and 
r9L6N256Wl, we verified that the steepness of the slope 
of the PDF over this range of overdensities is sensitive to 
the cosmological parameters employed. Note that the PDF 
does also not asymptote to Pv(A) oc A~'^, as predicted by 
Eq. IA5I For overdensities larger than the overdensity for the 
onset of star formation, the PDF is instead governed by the 
effective equation of state characteristic for the multi-phase 
star- forming gas. 

We fitted Eq. IA5I to our PDF over the range —1 ^ 
logj^g A ^ 2, constraining the parameters A and Co by the 
requirement that the total integral over the volume- and 
mass- weighted PDF must be normalized to unity. This yields 
values A = 3.038, 5o = 1.477, /3 = 3.380 and Co = -0.932. 
Note that the fit based on these parameter values still does 
not provide a good description of the PDF for overdensities 
login A > 2. 




Figure Al. Top panel: Volume- weighted PDF of the baryon 
overdensity A (per unit logjQ A) at z = 6 obtained from 
our default simulation r9L6N256 (black solid histogram). For 
comparison, the red solid line shows the four-parameter fit 
given by Eq. IA5I with the paramet er values taken from 
iMiralda-Escude. Haehnelt. fc ReesI l l2000t) . The blue dashed line 
shows our best fit of Eq. IA5I to the PDF. Bottom panel: Abso- 
lute value of the relative differences of the fits with respect to the 
PDF. In both panels, the vertical line indicates the overdensity 
corresponding to the onset of star formation. 
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